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Abstract 

The development and time evolution of a transport barrier in a magnetically confined 
plasma with non-monotonic, nonlinear dependence of the anomalous flux on mean gradi- 
ents is analyzed. Upon consideration of both the spatial inhomogeneity and the gradient 
nonlinear ity of the transport coefficient, we find that the transition develops as a bifurca- 
tion front with radially propagating discontinuity in local gradient. The spatial location of 
the transport barrier as a function of input flux is calculated. The analysis indicates that 

1/2 

for powers slightly above threshold, the barrier location Xb(t) ~ (D n t(P — P c )/P c ) , 
where P c is the local transition power threshold and D n is the neoclassical diffusivity . 
This result suggests a simple explanation of the high disruptivity observed in reversed 
shear plasmas. The basic conclusions of this theory are insensitive to the details of the 
local transport model. 
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I. Introduction 

Transport barrier physics is a central topic in ongoing research in magnetic fusion. 
By transport barrier, we refer to a region where anomalous transport is eliminated or 
very significantly reduced, so confinement is determined by neoclassical processes and 
macroscopic stability limits. It is important to note that a transport barrier need not 
necessarily be a thin layer (such the edge pedestal region in a standard High (H) mode 
plasma) , but rather can extend over a significant fraction of the plasma cross-section (as in 
the case of Enhanced Reversed Shear (ERS) modes) . As intimated above, the spectacular 
enhanced confinement characteristic of plasmas with transport barriers naturally presents 
a challenge to stability limits and the particle and exhaust control systems envisioned 
for advanced tokamak reactors. At the same time, the improved confinement is, of course, 
highly desirable. Hence, control of transport barriers (and transport in general) is emerging 
as an! important theme for present and future research in fusion plasma physics. One 
example of successful transport barrier control is the Ion Bernstein Wave (IBW) driven 
Core High (CH) mode achieved on the Modified Princeton Beta Experiment (PBX-M) 
tokamak. However, it seems fair to say that the science of transport barrier control is still 
in its infancy. 

Understanding the spatio-temporal dynamics of transport barrier is a necessary pre- 
requisite for engineering successful control techniques and algorithms. The many specific 
questions concerning the space-time characteristics of barrier evolution include: 

i. ) how can one predict the radial extent and location of a barrier? 

ii. ) what physics enters the criteria for barrier stationarity? 

iii. ) what is the rate of propagation of a non-stationary transport barrier? 

iv. ) what are the threshold conditions for barrier formation and relaxation? How much 

hysteresis is exhibited? 

v. ) how does a barrier respond to localized secondary external drive? For example, can 

localized deposition of IBW power be utilized to broaden an ERS transport barrier, 
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thus controlling profile peaking and reducing disruptivity? 
Here, we pursue several of these questions in the context of a simple, one-field model. 

Many physical problems 1 require the solution of the non-linear diffusion equation, 
where the diffusion coefficient D depends on the concentration of the diffused substance 
n. The power-low type of dependence on concentration: 

D~ n a 

leads to the formation of a propagating concentration front with the following behavior: 

n(x)^SKt)- 1 (l-(x/X(t)) 2 ) 1/a , if |*|<A(f); 
\0, if |x| > X(t), 

(see Fig.f). Here, X(t) = const ■ tV(2+a). j n this paper we consider a different type of 
nonlinear diffusion, namely one for which the diffusivity is a function of the gradient 
of concentration. Such a functional dependence occurs in magnetic plasma confinement 
devices, where the anomalous fluxes of high temperature plasma are driven by plasma 
microinstabilities with the characteristic length scale A ~ p t <C L r . Here pi is the ion 
Larmour radius and L r is the typical radial length scale of temperature (or density). 
In general, the fluxes of heat, particles and momentum are related to the temperature, 
density and velocity gradients, VT, Vn and VV, through the coefficients of the transport 
matrix. In the case of microturbulence driven transport, these coefficients are defined by 
the amplitude and correlation properties of turbulent fluctuations. The turbulence itself, 
however, is driven by the microinstabilities associated with the gradients of temperature 
and density, so the transport matrix is a function of VT, Vn and VV. Hence, the particle, 
thermal and momentum fluxes acquire a non-linear dependence on these gradients. 

The discovery of improved regimes of plasma confinement, such as H-mode, Very High 
(VH) mode and the improved core confinement modes 2 ' 3 ' 4 , strongly suggests, that this non- 
linear dependence is non-monotonic and the fluxes can actually decrease when the density 
or temperature gradients lie between certain values (i.e. negative differential diffusivity). 
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As a result, when the rate of fueling of plasma by particles or heat exceeds threshold value, 
a bifurcation to a new transport regime with higher values of Vn or VT takes place 5 ' 6 ' 7 ' 8 . 
The most plausible physical explanation of this phenomenon is based on the idea of the 
development of strong radial electric field shear, which stabilizes plasma instabilities and 
thus decreases the transport during the transition. In equilibrium, the shear of the radial 
electric field E r has contributions from the hydro dynamic plasma velocity V and pressure 
gradient dP/dr: 



where P; is the ion pressure and B is the magnetic field. The dynamics of perpendicular 
plasma velocity V are rather complicated 8 ' 9 ' 10 and thus are beyond the scope of this 
minimalist paper. Taking bulk velocity to be fixed, the force balance equation shown 
indicates that the changes in temperature and density gradients will be accompanied by 
changes in the shear of E r which, in turn, will change the anomalous transport. This 
constitutes a minimal model of a transport bifurcation. 

The purpose of this work is to describe the spatial and temporal development of trans- 
port bifurcation in the framework of a very simple, general one-field model by studing the 
geometry of flux surfaces in (x, <9 x n, V) space. The model involves the nonlinear diffusion 
equation with non-monotonic dependence of the local flux on the value of gradient. The 
model presumes a local transport processes yielding multiple gradients for a certain range 
of fluxes, only . In particular, no assumptions concerning the transport model or the bi- 
furcation mechanism are involved. We construct the solution in the form of propagating 
transport bifurcation front, which allows us to determine the final position of transport 
barrier from the condition of zero front velocity (i.e. stationarity). 

Although this work is discussed in the context of nonlinear transport phenomena in 
fusion devices, its results can also be relevant to other physical problems in which the 
diffusive fluxes depend on the gradients non-monotonically. Spinodal decomposition in 
alloys is one example of such a physical process 11 . Another one is turbulent mass and 
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heat transport in planetary atmospheres with zonal flows. In particular, analysis of the 
numerical simulations of Jovian atmosphere 12 show a suppression of radial thermal trans- 
port and generation of strong poloidal sheared flows when the radial drop in temperature 
(i.e. Rayleigh number) exceeds a certain critical value. This process is a classic example 
of a transport bifurcation, and is one which very likely determines the spatio-temporal 
patterns observed in the Jovian atmosphere. 

The dynamics of propagating transport barriers was already discussed in Ref.[13], 
where the spatially propagating front solutions describing the second order transition from 
an unstable transport regime to a stable one were introduced. The propagating bifurca- 
tion of our present model is different from that. It shares certain common features with 
front solutions in both the diffusion equation with the nonlinear diffusivity (D ~ n a ) and 
the bistable reaction-diffusion equation (i.e. Fitzhugh-Nagumo model) 13 ' 14 . The latter 
describes the dynamics of the first order transition from a metastble equilibrium state 
(or phase) to a stable one, corresponding to a global minimum of the Lyapunov function 
(effective potential energy) associated with the governing equation. The structure of our 
model equation is different from that classic example, but also has multiple stable equi- 
librium solutions in a certain range of input flux valu! es. As in the case of the reacti 
on-diffusion equation, the front solution in our system describes the transition between 
these equilibrium states (or transport "phases" ) . 

The remainder of this paper is organized as follows. In Section II we introduce basic 
equations and study the development of the transport barrier in the model with spatially 
dependent nonlinear diffusion. In Section III we study the space-time evolution of the 
propagating transport bifurcation and determine the stationary position of the transport 
barrier. In Section IV we discuss results, conclusions and their implications. 



II. Basic Model 

We consider the spatio-temporal dynamics of the simplest, one-dimensional transport 
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system which exhibits the property of negative differential diffusivity: 



d t n + d x T = Q(x), 



(1) 



where x (0 < x < a) corresponds to the radial coordinate in a tokamak and the flux V is 
given by the following expression: 



The detailed form of the function d x n) is not required here, it need only have the 
following very general properties 5 ' 15 : 
a. For any fixed value of d x n, $ decreases with x. The rate of this growth increases 
for x > x s corresponding to an increased level of anomalous transport at the plasma 



b. For any fixed value of x, is a function of d x n only, and has a characteristic S-curve 
shape (see Fig. 2) i.e. it increases for small values of <9 x n, goes through a maximum, 
decreases at intermediate values of d x n (corresponding to a negative differential dif- 
fusivity) and, finally, increases again for large values of the density gradient. The 
first interval of increasing $ is one of anomalous transport while the second stage of 
increase is determined by the collisional (neoclassical) diffusivity. Thus the slopes of 
$ (i.e. 5&/5\d x n\ ) in the two intervals of increasing flux are different. 
The contourlines of n') in the parametric plane {x, n'} are shown on Fig. 3, while the 
behavior of <E>, as a function of parameter n' only, is shown in Fig. 2 for several values of x. 

The linear differential operator L(d x ) accounts for smoothing effects, which are of 
higher order in the e ~ (A/L r ) 2 expansion. Thus L regularizes the small scale behavior of 
n. Here, A is the characteristic mixing length in the problem. It is defined either by the 
correlation length of the microturbulence or by the poloidal larmor radius (in the case of 
prevailing neoclassical transport). L is assumed to be an odd polynomial in d x , starting 
with the term d ■ d x . In the first approximation, though, the detailed structure of L is not 



T = 3>(x, d x n) — e L(d x ) * n. 



(2) 



edge; 
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crucial for describing front propagation, so long as this operator is dissipative. The source 
term Q(x) is a localized, even function of x, centered around x = 0. 
Let's analyze the stationary solutions of Eq.(l) with e = 0: 

d x ($(x,d x n))=0 (3) 

satisfying the following boundary conditions: 

n\ x=a = 0, r| x=0 = r . (3a) 

Eq.(3) can be easily integrated to yield: 

F = <S>(x,d x n). (4) 

This equation is easily analyzed graphically in the plane of parameters {x,n'}, where it 
defines a contourline <&(x,7i') = T = const. This contourline can be thought of as a plot 
of the function n' = n'(x,To) implicitly defined by $>(x, n') = T . It gives the solution of 
Eq.(3) with the boundary conditions (3a) in the following form: 

n (x) = / n'(s,T )ds. (5) 

J a 

Depending on the value of Fq, the function $ with the above-described properties allows 
the following three types of solution: 

a) . The solution u'^x^Tq) for small input flux values, r < T^, is shown in Fig. 4a. The 

parameter is defined as the local minimum of the S-shaped curve, which shows the 
dependence of the local flux on the density gradient at x = 0: 

F L = min$(0,n') for n < 0; 

b) . The solution n' 3 (a;,ro) for large values of flux, T > T Hl is shown in Fig. 4a. The 

parameter T# is defined as the local maximum of the flux curve at x = a: 

T H = max$(a,n') for n < 0. 
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In comparison with the previous case, this solution is characterized by a lower level 
of transport, which is given by the ratio r /n3(x,r ); 
c). The case with intermediate values of input flux, Th > r > T^, is shown in Fig. 4b. 
The equation n') = T can be inverted with respect to the parameter n' as n' = 
n'(x, To). This function, however, is multivalued in a certain range of x: < x cri (To) < 
x < x cr2 (r ) < a because of the S-shaped form of the corresponding contourline. As 
a result, we can identify three branches of n'(x, T ), which will be further referred to 
as n'^x), i = 1,2, 3, satisfying > n[(x, To) > n 2 (x, To) > n' 3 (x, To). 
In case a) or in case b), the solution of the reduced system (3), (3a) given by Eq.(5) ap- 
proximates the stationary solution of Eq.(l) (when the source term is absent) with the 
relative accuracy e everywhere in x except for narrow boundary layers at x = and 
x = a. In case c), however, this solution is not well defined because of the presence of 
several branches of n'(x,T ). Nevertheless, it is possible to build a composite solution 
n c (x,To) = f n' c (s,To)ds, where the function n' c consists of two or more smooth parts 
from the different branches of the multivalued function n'(x,T ) which are separated by 
jump discontinuities. In this solution, the branch with larger value of \n' c \ corresponds to 
the transport barrier. Such a composite solution is a good approximation of the stationary 
solution of Eq.(l) everywhere in x, except for the narrow boundary layers in the vicinity 
of the jump discontinuities and at the bou! ndaries x — > 0, a (see Fig. 5). These provide 
smooth transition between the branches corresponding to different transport regimes. The 
locations of these boundary layers, as well as the possibility of their propagation as dynam- 
ical fronts will be discussed in the next Section. One should observe, that the composite 
solution can be constructed only out of the first and third branches of the multivalued 
function n'(x,T ). Its second branch n' 2 (x,r ) belongs to the area of {x,n'} plane which 

is characterized by negative differential diffusivity: 

d$(x,ri) 



$'(z) 



< 0. 

n'=n' 2 



dn' 

As a result, it is unstable with respect to small perturbations n(x, t) = n(x, t) — riiix, To). 



The linearization of Eq.(l) yields 



d t n = d x ( d x n ) + eL(9 :E ) * n. 



The spatial scale of the perturbation, /, is assumed to be small i.e. I <gi L r , but still larger 
than \feL r . The above equation can be written in the following form: 



For < 0, it produces unstable solutions with small scales growing fastest, in the absence 
of regularization by L. In particular, for an S-shaped $>(x, d x n) the instability related to 
<&' < value of d x n is what drives the propagation of the bifurcation. The analogous 
drive for a super-critical bifurcation front is the instability familiar from the Fisher and 
time dependent Ginzburg-Landau equations, discussed (in the context of transport barrier 
dynamics) in Ref.[15]. The width of a transitional boundary layer corresponding to a 
jump in derivative of the composite solution can be estimated from the fact that this layer 
matches the n' x and n' s branches of n'(x, Fq) through the values of d x n corresponding to a 
negative differential diffusivity: (d$>/dn')\ n >=d xn . This is possible only when the change of 
the nonlinear flux function 5<E> accross this layer is balanced by the d! issipative operator 



This condition immediately yields the following estimate of the boundary layer width 



The structure of the stationary solutions of Eq.(l) is such that for continuously chang- 
ing r , transitions between different transport regimes occur in the form of bifurcations. 
As the value of r is increased above Tl, the solution ni(x,F ) with r < T^, continu- 
ously evolves into the lower branch ni(x,r ) with r > Fl- This solution continues to 



d t h « $' din. 



eL: 



5$ ~ r ~ e L(d x ) *n ~ ed 



d x n 
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change smoothly with the further increase of T until the flux reaches the threshold value 
corresponding to the local maximum of the flux curve at x = 0: 

Tf r = max $(0, ri) for ri < 0. 

When To > T^ r , the function ni(x,To) doesn't exist for < x < x cri (To) < a, so a new 
solution should be sought either in the form of a composite solution n c (x,T ) or in the 
form of the branch n 3 (x,r ). Hence, for r ~ r^"", a small increase of flux results in a 
significant jump of — d x n(x,To) to a much higher value in the interval < x < x cri (To). 
This bifurcation can be described as a formation of transport barrier. The width of this 
barrier will be set by the stationary position of the transitional boundary layer connecting 
the zones of enhanced and suppressed transport. Apparently, this width can't be less than 
x cri (T ) and more than x cr2 (T ). In principle, when a; cr2 (r ) > a it can cover the whole 
range of x. The issue of the barrier width will be addressed in the Section III in more 
detail. For the values of flux decreasing from some starting value To exceeding Th, the 
profile of n(x, T ) exhibits a similar bifurcation of the solution n 3 (x, r ) at T = T^ r . This 
quantity is defined as the local minimum of the flux curve at x = a: 

T^ r = mm ^(a,ri) for ri < 0. 

In this case, a transition to the regime with high transport at x cr2 (To) < x < a takes 
place. 

It is clear from the above that the presence of two locally stable branches of the solution 
n(x,T ) for Tl < r < Th results in hysteresis. For example, when r is increased, the 
transition from ni(x, T ) to a profile with a lower level of transport occurs at r = < 
Th, while the solution with the profile n^{x, Tq) will bifurcate to a profile with a higher level 
of transport at r = T l ^ r >T L , when the flux is decreased. As a result, when T l ^ r > 
is satisfied, the profile will return to the mode with high level of transport ni(x,T ) at 
a value of the flux lower than the one that is required for the bifurcation from ni(x, Tq), 

10 



when the flux is increased. Similar hysteresis behavior is also possible for < . In 
this case, the bifurcation to the profile with a lower level of transport will occur at the flux 
value r^ r , when the left branching point crosses the left boundary: ^^(To) = 0. For the 
interm ediate values of flux > r > T^ 17 ", the position of the transitional boundary 
layer Xb(To) is somewhere between x cri (To) and £ cr2 (To). When the flux is decreased 
below T^ r , the transition to the regime with a high level of transport will occur at the 
moment when Xb(To) crosses the left boundary: Xb(To) = > x cri (To). Apparently, this 
will happen at the value of flux which is lower than r^ r . These bifurcation scenarios are 
illustrated in Fig.6a,b. In a simplified model with two different linear diffusivities D an 
and D neo for \n'\ < \n' crit \ and \n'\ > \n' crit \ respectively, the hysteresis ratio (the ratio of 
the input power necessary for the transition to a higher confinement mode to the power 
at which the transition back occurs) scales as the ratio of the anomalous to neoclassical 
diffusivities i.e. Y l ^ h /^h\^l ~ D an /D neo . 

III. Spatial dynamics of transport barrier. 

When analyzing the dynamics of the solutions with transitional boundary layers, we 
can neglect the spatial dependence of the nonlinear flux function $ as long as the tran- 
sitional layer width remains small i.e. A x ~ e 1 / 2 L r <C L r . Let's consider the following 
modification of Eqs. (1),(2) : 

d t n + d x T = 0, where V = $(d x n) - edd x n. (6) 

The function $(n'), which is shown in Fig. 7, depends on the variable n' in the same way 
as does the function n') with a fixed value of x. For function $ independent of x, the 
threshold values of r defined in Section II satisfy the following relations: 

rf r = r H , and r% r = r L . 



The coefficient d in the last term of the Right Hand Side (RHS) of Eq.(6) is of the order 
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of (Ll<fr)/d x n. Equation (6) has the following boundary conditions: 

<f>{d x n\ ) = T , n\ a = 0, d 3 x n\ {0 , a} = 0. (7) 

We are interested in the values of flux which allow for multiple stationary solutions of 
Eq.(6): Th > T > T^. For these values of r , there are two stable, stationary solutions: 

ni(x, r ) = n'i( r o) (x - a), and n s (x, F ) = n' 3 (T Q ) (x - a), 

where n' 13 (To) are the first and the third roots of the equation $(n') = To (see Fig. 7). If 
one introduces the parameters n' L H corresponding to the positions of the local minimum 
and maximum of the function $(n'): <&{n' L H ) = T^ ^, then the inequality —n' 3 (T ) > 
—n' L > —n' H > —u'^Tq) > is valid. We look for an asymptotic solution of Eq.(6) with 
the boundary conditions (7) corresponding to the limit e — > 0. This solution should consist 
of two pieces separated by a boundary layer at the point Xb- Each of these pieces lies on 
a separate, stable branch of the S'-shaped flux curve $(n'): 

J —d x n(x, t) > —n' L , for < x < Xb — A, 
[ —d x n(x, t) < —n' H , for Xb + A < x < a. 

Here, A is a free parameter which is small compared to L r , but exceeds the width of the 
layer: A > A b . Everywhere outside of the transitional layer: Xb — A < x < Xb + A, the 
solution is assumed to be a slowly varying function of x with the characteristic length scale 
L r , comparable to the system size a. 

In a stationary state, the necessary condition for the existence of this solution can be 
obtained from the integration of the expression for the flux conservation r = &(d x n) — 
edd^n multiplied by d^n. The integration yields: 

"d x n\ c 
I d x n\, 

When e — > 0, the above relation has the following geometrical interpretation: the total 
area between the graphs of V = r and V = $(n') over the interval between the points of 
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-r d x n 



rd x n\ a a 

/ ^{n')dn' + -{dlnf (8) 



their intersection (ji'^Tq), n' 3 (T )} is equal to zero, (see Fig. 7). For every curve $(n'), this 
relation specifies a single value of r . It will be hereafter referred to as Y M - The condition 
yielding Ym is, of course, closely related to the Maxwell construction criterion for phase 
equilibrium. 

For r 7^ Tm, a stationary composite profile does not exist . Nevertheless, in this 
case a time dependent composite profile can be found. It consists of two parts, one with 
— d x n < —n' H and the second — d x n > —n' L , separated by a propagating boundary layer at 
x 7a Xb{t) (see Fig. 8). In order to construct such a profile, let's consider a model nonlinear 
flux function: 

( d a n' ', - anomalous transport for —n' < — n' H — 5; 
3>(n') = < d n n\ - neoclassical transport for —n' > —n' L + 5; (9) 
y g(n'), -suppressed anomalous transport for — n' L — 5 < —n' < —n' L + S. 

In this very simple model of nonlinear flux behavior in a tokamak plasmas, the linear 
diffusion coefficients d n and d a correspond to the low level, neoclassical transport in the 
core (0 < x < Xb) and anomalous transport at the edge (xb < x < a). The function g(n') 
smoothly connects the linear branches of $(n'). Its local maximum and minimum are at 
the points n' H and n' L , correspondingly. In addition, the inequality Q(n' H + 5) > Ym > 
&(n' L — S) is assumed to be satisfied, where Ym is the value of flux obtained from the 
Maxwell construction for $(n'). For small e, the solution with a transitional layer has the 
following form: 

a _ / d n din, and — d x n < — n' L — 5 < 0, for < x < Xb(t) — A; , , 
\d a d1n, and > — d x n > — n' H + 5, for a > x > xj,(t) + A, 

with the following boundary conditions: 

$(a x n)| =r , n\ a =o. 

The first two matching conditions at the boundary layer are rather straightforward i.e. 
the continuity of n: 

™L-A=< + A+°( A )> (11) 
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and the continuity of flux: 

F L-A= d nd x n\ Xb _ A = r| Xb+A = d a d x n\ Xb+A +0(A). (12) 

An additional matching condition is obtained by integrating Eq.(6) across the transitional 
layer with the weight d x n, 



Td x n 



x b +A rd x n 



/ $(n')dn' + 0(A). (13) 

J d x n\ x a 



x b -A 

As noted above, this relation is an example of a Maxwell construction 16 ' 14 ' 18 , and is related 
to the condition for the coexistence of two transport regimes ("phases"). Its geometrical 
interpretation is identical to that of Eq.(8), shown in Fig. 7. 

As A, e — > 0, the above system results in the following linear problem with the surface 
of discontinuity at x = xi(t): 



(14) 



j d t n = d n d^n, for < x < Xb(t), 
\ d t n — d a din, for a > x > Xb(t), 

n\ a = 0, $(d x n\ ) = F , 
n\x b -o = n\x b +o, d n d x n\ Xb - = d a d x n\ Xb+0 = —T Xb , (15) 

x b +0 fd x n\ Xb+0 

T Xb d x n = / $(n')dn'. 

*b-0 Jd x n\ Xb - 

Note, that flux continuity in combination with the last condition of transport "phase" equi- 
librium are equivalent to the definition of the quantity r^, so the flux at the transitional 
layer is fixed at T Xb = Fm- As a result, these two matching conditions can be rewritten in 
the following form: 

— d n d x n\ Xb -o = T M , —d n d x n\ Xb+ o = r M - (15a) 

When the flux T on the left boundary coincides with T M , the system (14), (15) gives 
a trivial solution for a^: 

Xb = const, < Xb < a. 
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Let's now consider the situation when T slightly exceeds T M - In that case, Xb(t) can't 
be constant. Otherwise, the asymptotic (in time) solution of Eq.(14) for x > Xb would 
be: n(x) = — (ITm/c^) (x — a), giving the value of n at the transitional layer n(xb) = 
(Tm /d a ) (a — Xb) = const. However, this value cannot be matched with the t — > oo asymp- 
totic of the solution for < x < Xb, which increases with time as n « (Tq — TM)t/xb- 
In principle, the system (14), (15) can be solved exactly. Here, we seek an approximate 
solution, which describes a slow propagation of the transitional layer, i.e. which satisfies 

d t x\ < d n . 

In practice, this is equivalent to requiring that the barrier propagation velocity Vb = dtXb 
satisfy Vb < d n /xb- The diffusivity d a exceeds d n so the relaxation time r a ~ a 2 /d a for 
the solution at the edge (x > Xbit)) is small compared to r n ~ a 2 /d n for the solution in 
the core (x < Xbit)) to develop. Hence, the solution for x > Xb(t) can be taken to be 
stationary, i.e. 

n(x) = (T M /d a ) (a - x), for x > x b (t). (16) 
For x < Xbit), we can make the following substitution: 

( \ r ° i ( r o - Tm) x 2 
n(x) = 'T n X+ 2xW) + /(X ' ^ (17) 

where / is a new unknown function satisfying: 

f) f -j a2 f (r -r M ) , (Tq-Tm) x 2 d t x b (t) 

d t f = d n dj + xb{t) + — (18) 

d x f\o = d x f\ Xb = 0. 

Here we use the boundary conditions defining the flux at the left boundary, x = 0, and 
at the transitional layer x = Xb, only. According to the assumption of the slow front 
propagation, the last term in the RHS of Eq.(18) can be neglected. As a result, we obtain 
the following approximate solution for x < Xb(t): 

< * r « , ( r o-r M ) x 2 r dt 
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Matching of this solution with the solution for x > Xb yields the equation for x\,{t): 



where C is a constant of integration. The assumption of slow front propagation is satis- 
fied for values of flux T close enough to Tm- (^o — ^m)/^m <C 1. The profile n(x,t), 
corresponding to the above solution is shown in Fig. 9. This simple analysis allows us to 
conclude, that when the transport barrier is created in a certain range of x, it will spread 
as long as the condition T > Tm is satisfied. For d a ~^> d n , the rate of its propagation will 
be mainly defined by the low diffusion rate in the core: Xb(t) ~ \J ((T — Tm)/(2Tm)) d n t. 
When the x-dependence of the nonlinear flux function is taken into account, our analysis 
may be considered as quasi-local. The conclusions should not significantly change as long 
as the transition layer width remains smaller than L r . In the opposite limit of d t x\ S> d n , 
no front solution exists. 

The description of the bifurcation scenarios started in Section II can be completed 
now. For example, when > r > r^ r , the solution has multiple branches at < 
^cri(ro) < x < x cr2 (ro) < a, the ambiguity in the radial extent of the low transport 
zone can be resolved. Its boundary is at the point x, where the local front velocity d t Xb, 
obtained from r and the local value of the equilibrium flux rjif(x), is zero. This takes 
place, when the local "phase" equilibrium condition is satisfied: To = Ym{x). Apparently, 
TM(x cri ) < r and TM(x cr2 ) > r , so such a point can be found in the interval (x cr2 , x cr2 ). 

IV. Summary and Conclusions. 

Experimental results and theoretical models of the anomalous heat and particle trans- 
port in magnetically confined plasmas suggest that they possess the important property 




(20) 



This can be easily solved, giving: 




(21) 
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of negative differential diffusivity for a range of temperature and density gradients. As a 
result, when the fueling rate of the plasma exceeds a certain threshold, the profiles undergo 
a bifurcation to a confinement regime with much steeper gradients (i.e. transport barrier 
forms). In this paper we have demonstrated a simple dynamical model of such barrier 
dynamics. The main conclusions of this paper are summarized below. 

i). the spatial localization of the transport barrier is determined by the structure 
of the nonlinear flux function in {x, d x n} parameter plane. The shape of this function 
allows for multiple solutions for the profile n(x). The resulting profile is represented by a 
single solution branch or by a combination of branches corresponding to neoclassical and 
anomalous transport connected through transition layers. The stationary position x\> of 
these layers is given by an argument similar to that used in the Maxwell construction. 
Specifically, the layer is stationary at the point where the total flux through the system 
coincides with the local value of the Maxwell flux: 



ii). the spatial dynamics of the transition evolves in the form of a propagating tran- 
sition layers. For d t x\ < d n , the speed of their propagation is defined by the neoclassical 
diffusivity and the difference between the input flux T and the local Maxwell flux: 



It seems likely, that for values of flux r significantly different from Fm, |To — ^m\ ~ Tm, 
the front velocity is still restricted by the slowest (neoclassical) diffusion rate. This follows 
from the observation, that for < x < Xb profile evolves on a time scale x\jd n , which 
should be comparable with the time scale of the front propagation, i.e. Xb/dtXb- 



T = T M (x). 



The flux Tm(x) is given by the following equation: 




where ${n'< 
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iii). the transition between different transport phases exhibits hysteresis. Even in the 
simple model considered in our paper, several types of hysteresis are possible. It occurs 
both when the flux To is increased beyond r^ r and subsequently reduced, and when To 
is increased above Th > r^ 17 " and later reduced. 

The transport barrier propagation phenomenon described here is somewhat similar to 
the spatial evolution of a first order transition from a metastable phase to a stable phase in 
a non-equilibrium thermodynamic system. The condensation of overheated vapor 17 or the 
spinodal decomposition 12 are examples of such processes. Our results may also be relevant 
to the anomalous transport in convective systems, such as planetary atmospheres, where 
the turbulence coexists with large scale, self-consistently generated zonal flows. 

As mentioned in the introduction, understanding the space-time response of transport 
barriers is critical to transport control, which is a likely element of any realization of an 
advanced tokamak. In this regard, one particularly compelling motivation for transport 
control is to reduce and minimize the disruptivity of core transport barrier plasmas thus 
allowing full exploitation of their improved confinement. The disruptions which occur in 
core barrier plasmas are almost certainly a consequence of the dramatic steepening of the 
core pressure gradient due to barrier formation at high heating power. In the context of 
the simple model discussed here and in Ref.[15], the core pressure gradient is determined 
by: 

a) , the source strength and (fixed) anomalous diffusivity, 

b) . the transition layer location Xb(t), which supplies the "boundary" condition. 

In simple terms, our model suggests that the global profile evolution will resemble swinging 
doors attached to a sliding hinge, located at Xbif). Thus, the energy content can either be 
increased by Xb(t) moving outward or by increasing P' core , eventually leading to violation 
of macroscopic stability limits. Hence, the dynamics and stationary value of Xb(t)p\a,j an 
important role in determining disruptivity! One of the central results of this analysis is 

1 /2 

that the hinge cannot slide very fast (i.e. Xb(t) ~ (d n t(P — P c )/P c ) ). This situation 
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is exacerbated by the radial dependence of P c , which increases sharply with r as s goes 
positive 15 . Since the barrier cannot expand quick enough, P' core surges due to the combined 
effects of strong fueling and low transport, which in turn steepens the profile and leads 
to disruption. A straightforward corollary of this hypothesis is that disruptivity should 
decrease for broader barriers. This is borne out in Doublet III-D (DIII-D) Negative Core 
Shear (NCS) discharges, where an L— >H transition is triggered (thus broadening the de- 
position profile), and Weak Negative Shear (WNS) plasmas, in which the q(r) profile is 
rather flat, thus reducing the radial gradient in P core (r). It should be noted, however, that 
control via an H-mode edge is undesirable, since the beneficial properties of the L-mode 
are lost. Alternative control schemes, which exploit Radio Frequency (R! F)-driven shear 
layers to broade n the transition barrier are being examined 18 and will be discussed in 
future publications. 
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